New method to study stochastic growth equations: a cellular automata perspective 
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We introduce a new method based on cellular automata dynamics to study stochastic growth equa- 
tions. The method defines an interface growth process which depends on height differences between 
neighbors. The growth rule assigns a probability Pi(t) = p exp[« Ti(t)] for a site i to receive one 
particle at a time t and all the sites are updated simultaneously. Here p and k are two parameters 
and is a function which depends on height of the site i and its neighbors. Its functional form 

is specified through discretization of the deterministic part of the growth equation associated to 
a given deposition process. In particular, we apply this method to study two linear equations - 
the Edwards-Wilkinson (EW) equation and the Mullins-Herring (MH) equation - and a non-linear 
one - the Kardar-Parisi-Zhang (KPZ) equation. Through simulations and statistical analysis of the 
height distributions of the profiles, we recover the values for roughening exponents, which confirm 
that the processes generated by the method are indeed in the universality classes of the original 
growth equations. In addition, a crossover from Random Deposition to the associated correlated 
regime is observed when the parameter k is varied. 
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INTRODUCTION 



roughness u>(L,t), which is defined at time t as 



Discrete computational growth models have been 
largely investigated along the last decades, due to the 
great interest in describing various features of interface 
growth phenomena, observed in a wide range of phys- 
ical processes jH 0. Q. As examples of kinetic rough- 
ening models we can mention the Eden model Q|, the 
ballistic deposition model (BD) Q and some solid-on- 
solid growth models in which correlation mechanisms are 
present, such as surface relaxation height difference 
restriction curvature restriction |8l and surface diffu- 
sion [9, 10]. From an experimental point of view, growing 
interfaces can be generated, for example, by Molecular 
Beam Epitaxy (MBE) and vapor deposition over cold 
substrates (see Refs. 0,01 and references therein). 

The methodology used to study interface growth phe- 
nomena have also been applied to investigate Cellular 
Automata (CA) [ijj, an immense class of computational 
models which describe many phenomena in a wide sort 
of scientific subjects. By means of an accumulation 
method 0], usually one can map the evolution of a CA 
into a growing profile and then apply the tools used to an- 
alyze such systems. This method has been used recently 
to study deterministic an d probabilistic CA 

In the analysis of interface growth one is generally 
concerned about the temporal behavior of the interface 
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where L is the size of the substrate and h(L,t) is the 
mean height of the generated profile. 

For growing systems it is known that the roughness 
grows initially with time f as a power law, defining the 
growth exponent, (3. After a saturation time t x , the in- 
terface reaches a stationary regime and the roughness 
saturates. Both the saturation roughness, uj sat , and sat- 
uration time, t x , depend on the system size, L, as a power 
law, defining the roughness exponent, a, and the dynamic 
exponent, z, respectively. The roughness of the surface 
follows the Family- Vicsek scaling law 0] 



Lu(L,t)~L a f 



(2) 



where the scaling function f(u) behaves as f(u) ~ vP for 
u <C 1, and f(u) = constant for m>1. It follows that 
f3 = z/a. 

A set of values for the roughening exponents, in a given 
dimension, specifies a Universality Class (UC). Thus, if 
two or more processes have the same exponents values, 
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one can say that they belong to the same UC, which 
means that their underlying dynamics have the same 
symmetries and conservation laws. 

Based on these symmetries, one can write a stochastic 
differential equation and associate it to the given UC. 
This continuum approach often provides an analytical 
treatment for growing interfaces and exact values to the 
roughening exponents can, sometimes, be obtained. 

In this paper we consider two linear equations and 
one non-linear one. They are, respectively, the Edwards- 
Wilkinson (EW) equation 



tA7 2 /i(x,i) +rj(x,t) , 



dh(x, t) 
Ft 



the Mullins-Herring (MH) equation [1J, [l 



(3) 



dh(x, t) 

at 



-KX7 4 h{x,t)+r}(x,t) 



(4) 



and the Kardar-Parisi-Zhang (KPZ) equation [19 ]. 



9 -^- = ,V 2 Mx, t) + \ (V/i) 2 + t) . (5) 

Here, /i(x, t) is the local height of the profile, assumed to 
be continuous (eventual hangovers are to be ignored) , and 
7?(x, t) is a gaussian noise, with zero mean and correlation 
given by 



(r? (x, t) V (x', t')) = 2DS d (x - x') S(t- t') . (6) 

The linear equations can be solved exactly and the val- 
ues of the roughening exponents can be determined. For 
the UC associated to the EW equation CQ| , the exponents 
for a <i-dimensional substrate are 



and 



(7) 



This equation can be associated to the Random Depo- 
sition with Surface Relaxation model (RDSR) ;0|, where 
the particles, after being deposited in a random position 
on the lattice, are allowed to relax to the local minimum, 
considering the nearest neighborhood of the chosen site, 
and sticks irreversibly to the aggregate. 

The solution for the MH equation © yields U03 
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On the other hand, non-linear equations have no gen- 
eral solutions for the roughening exponents. Neverthe- 
less, for the KPZ equation d^fc with d = 1, renormaliza- 
tion group theory provides [19] 
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(9) 



This equation can be associated to two discrete mod- 
els: the Ballistic Deposition (BD) |j| and the Restricted 
Solid-on-Solid growth model (RSOS) 0. In the BD, the 
falling particle sticks to a vertical position such that 



hi(t + l) = 



fti(t) + 1, - (10) 



In the RSOS model, at each time step, a random po- 
sition is chosen in the lattice and its height is increased 
by one unit, provided the restriction \Ah\ < m is obeyed, 
where Ah = hi — hi±\. 

Crossovers between distinct UCs is a topic of great 
interest as well. Much attention have been given to 
competitive growth models, where two different types 
of particles are deposited, one with probability P and 
the other one with probability (1 — P), allowing several 
combinations to crossover between different growth pro- 
cesses (2lll22l|. In particular, models where some kind of 
correlated growth occurs with probability P and Random 
Deposition (RD) with probability (1 — P) have been in- 
vestigated recently, by means of simulations and scaling 
arguments [23L l2~3 | . A crossover from random to corre- 
lated growth in the RSOS model have also been studied 
recently by Aarao Reis |24[ and a crossover from KPZ to 
EW class has been obtained by da Silva and Moreira in 
a growth model with restricted surface relaxation [25] . 

Our purpose is to develop a method [2f| to study those 
equations, based on CA dynamics. By using a simple dis- 
cretization scheme, we obtain numerical solutions for the 
roughening exponents without actually having to solve 
it, which means that the method can be applied to equa- 
tions whose solutions are not even known. In section 
HTl we outline the main features of the method such as 
definition of the parameters and their range of interest, 
discretization schemes and a brief description of the algo- 
rithm. In section ITTTI we show numerical results obtained 
for the roughening exponents corresponding to each UC 
considered. In section HVl we present a description of the 
crossover from Random Deposition (RD) to correlated 
growth obtained by variation of K. Finally, in section Fvl 
we draw some conclusions and perspectives for further 
works. 



and z = 4 
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This equation is associated to a deposition model known 
as Growth with Surface Diffusion (GSD), where the 
newly arrived particle seeks the position in the neigh- 
borhood where the bonding is maximum [flllicj. 



II. THE METHOD 

Consider a one-dimensional lattice of size L, initially 
flat and with periodic boundary conditions. In each time 
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step, all the sites are simultaneously visited and site i 
receives a particle with probability Pi(t), given by 

Pi (t) = p e KT ^ . (11) 

Here < p < 1 and k > are two parameters, fixed 
throughout the evolution of the interface. The former 
is related to the growth speed and the later can be as- 
sociated to an inverse of temperature if one considers 
an analogy with vapor deposition: the synchronous up- 
date scheme can be thought as attempts for deposition 
of vapor particles on a cold substrate, allowing the de- 
velopment of growing structures. Ti(t) is the kernel, a 
function which depends on heights of site i and its neigh- 
borhood, at time t. Its explicit form will be given by 
the discretization of the deterministic part of the growth 
equation we are intended to study. 

In the case of the EW equation Q , the kernel is given 
by the discretization of the Laplacian V 2 /i, 

r*(i) = h l+1 (t) + hi-!® - 2h.it) . (12) 

For the MH equation the kernel follows from the 
discretization of the negative of the fourth spatial deriva- 
tive, -V 4 /i, 

Ti(t) = -[6hi{t)+h i+ 2{t)+hi-2{t)] + 

+ 4[ft i+ i(t) + fc i _i(t)] . (13) 

Finally, the discretization of the square of the gradient, 
(Vft.) 2 , and the Laplacian yields the kernel for the KPZ 
equation (§J, which is 

Ti(t) = -[h i+1 (t)-h i -i(t)] a + 

£ 

+ [h i+ !(t)+hi- 1 (t)-2h i (t)] , (14) 

where e > is the parameter which controls the non- 
linearity strength: large e implies a small contribution of 
the non-linear term, and conversely. 

So far we have explicitly considered only the deter- 
ministic part of the growth differential equations <EJ| to 
©. Nevertheless, the stochastic nature of such equations 
contained in the gaussian noise r?(x, t) is simulated in our 
method by the probabilistic character of the growth pro- 
cess. We have done a rigorous study of the symmetry 
and decay properties of the height distributions, which 
corresponded to those of a gaussian distribution. 

Note that, in the way we have defined the method, one 
can eventually obtain Pi(t) > 1. For this situation, we 
impose the condition 

Pi(i) > 1 => Pi {t) = 1 => hi(t+l) = hi{t)+l . (15) 



Hence, given a pair of values (p, k), there is a maximum 
kernel value, T max (p, k), for which 

> T max Pi (t) = 1 . (16) 

Making pi (t) = 1 in equation ifTTJl and having in mind 
the fact that Ti(t) is an integer by definition, we find 



r ma x(|0, «) = int ( - - lnpj . (17) 

Basically, for each time step t the algorithm does the 
following: (i) calculate r,*(t); (ii) ask whether Ti(t) > 
F m ax', (iii) if so, a particle should be deposited in site 
i; if not, a random number r is taken in the range [0,1) 
and Pi(t) is calculated: if r < Pi(t), a particle should 
be deposited on site i; (iv) repeat steps (i) to (iii) for 
i = 1, . . . , L; (v) those sites that should receive a particle, 
have their heights simultaneously increased by one unit. 

For consistency between the algorithm and the def- 
inition of the method, equation 111 111 , we must have 
^max > 1) because T max — 0, according to the algorithm, 
would make the interface grow flat since the beginning of 
the process and no scaling features would be observed. 
This condition restricts the range of values of K that are 
to be considered, 

r max > 1 K < In . (18) 

We fixed p —\ throughout our simulations so that, in 
a flat surface, each site has probability p — \ to receive 
a particle. Further analysis of the displacement rate of 
the mean height showed that the parameter p is related 
to the growth speed of the interface. 

III. ROUGHENING EXPONENTS 

In this section we show the results obtained for the 
roughening exponents, for each one of the three UCs. In 
a lattice of size L = 10 4 , we initially fixed K = 1CP 1 and 
averaged the results over 50 independent samples. 

The results obtained for the roughness of the interfaces 
generated by the application of the method to the EW 
and MH equations are shown in figure Note the good 
agreement between the growth exponent values obtained 
in our simulations and the predicted ones, equations JJJ 
and JHJ with d= 1. 

In the application of the method to the non-linear case 
of the KPZ equation, our simulations have shown that 
in the situation where the contribution from the Lapla- 
cian vanishes (e — > 0) and the system evolves following 
a pure non-linear dynamics, one obtains persistent un- 
evennesses growing between "valleys" and "hills", because 
the square of the gradient induces the development of 



4 




Figure 1: Log-log plots of the roughness as function of time, 
for p — 1/2, k = IO -1 , L — IO 4 and averaged over 50 inde- 
pendent samples. The circles represent the data for the EW 
class ((3 = 1/4 expected), while the triangles are the data for 
the MH class ((3 = 3/8 expected). Least squares method was 
used in order to determine the indicated numerical values for 
fl- 



ench structures. We can see in the top frame of figure 
|2l the evolution of a profile generated in the pure non- 
linear case, where the set of probabilities reaches a stable 
configuration such that Pi — constant for all i: sites with 
larger heights have pi = 1 and sites in the valleys have 
Pi = p (these valleys are symmetric, i.e. both sides have 
same slope). Thus, when that stage is attained by the 
system, the roughness, which is a measure of the width 
of the interface and can be thought as the difference be- 
tween maximum and minimum heights, grows linearly 
with time and we get (3=1. This behavior is shown in 
the bottom of figure It is worthy to mention that this 
behavior is independent of the size of the system and 
occurs for all n, the crossover from (3 = 1/2 to /3 = 1 
happening for larger times as n decreases. 

In the top of figure we show the profile generated 
when one sets e with small values (e = 2), which corre- 
sponds to a small, but non- vanishing, contribution from 
the Laplacian term. In this case the method produces 
crystallized patterns with asymmetric hills and valleys, 
in which all sites i have Pi(t) = 1 for all t (after a cer- 
tain transient time) and thus the roughness no longer 
changes. In the bottom of figure we show the behavior 
of the roughness for several system sizes and, as one can 
see, for larger systems the roughness grows initially with 
(3 w 1/3 before the frozen configuration is reached, while 
smaller systems can saturate before that. We believe 
that these structures are attractors among the possible 
configurations of heights in the profile. So, for small e, 
the system always reaches such absorbing configurations. 
In our simulations the size of the system is restricted to 
L < 10 5 so that we must make e > 4 in order to have 
the standard behavior of the roughness and to avoid such 
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Figure 2: In the top frame we show the evolution of the profiles 
by changing the color of the particles each 10 2 time steps, 
for 10 2 < t < 10 3 , in the limiting case of e — > 0. In the 
bottom we display the log-log plot of the time behavior of the 
roughness for L = 10 3 and several values of k, between 10 _1 
and 10~ 4 . Note that for smaller values of k the crossover 
(3 = 1/2— >/3 = l occurs at larger times. 



anomalous configurations. 

We show in figure 0}a the results obtained for the 
roughness, with e = 5. Again, as one can see, a good 
agreement with prediction was obtained for the growth 
exponent (3. By making e larger, a crossover from 
(3 = 1/4 to (3 = 1/3 regimes is observed, as one can 
see in figure 0}b, where we have e = 10 for a system of 
size L = 10 5 . This crossover occurs because in the be- 
ginning of the growth process, when the roughness is not 
large enough, the quadratic term of the KPZ equation is 
much smaller than the Laplacian, which dominates and 
provides (3 « 1/4, As the roughness increases, the non- 
linear term becomes dominant and we get (3 ~ 1/3. Of 
course, if we let e — > oo, this crossover does not occur 
anymore and we recover the results obtained for the EW 
class. This crossover has been previously obtained by 
da Silva and Moreira in a deposition model with re- 
stricted surface relaxation, where particles can relax only 
within a given distance s. If a minimum cannot be found 
in this range, the particle evaporates in a way similar to 
the RSOS growth model fjj. An s-dependent crossover 
from (3 = 1/4 to (3 = 1/3 was obtained by the authors. 
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Figure 3: In the top frame we show the interface evolution for 
the case of e = 2, for 10 6 <t< 1.5xl0 7 , where we change the 
color of the particles each 10 6 steps. In this case the system 
is very close to a crystallized pattern. In the bottom we have 
k, = 1CP 1 and several values of L in the range 25 < L < 10 3 .- 
larger systems reach crystallized patterns while smaller ones 
can saturate before that happens. 




Figure 4: Log-log plots of the roughness as function of time, 
in the application to the KPZ equation, where (i = 1/3 is 
expected. In the top, we have e = 5 for a system of size 
L — 10 4 . In the bottom, e = 10 and L — 10' , where a 
crossover from f3 = 1/4 to [3 = 1/3 is observed. We have 
drawn the functions it) ~ t 1 ' 4 (dashed line) and u ~ t 1 ' 3 
(dotted line) for comparison. 



IV. CROSSOVER FROM RANDOM TO 
CORRELATED GROWTH 

As we vary the parameter n, making it smaller, we 
identify a crossover from RD regime ({3 — 1/2) to the 
corresponding correlated process. It was found that this 
crossover does not depend on the system size L. 

Defining the crossover time as t c , we see that t c , iv sa t 
and i x are all functions of the parameter k, in a power 
law fashion: 



To determine the roughness exponent a and the dy- 
namic exponent z, in order to have a complete charac- 
terization of each UC, we varied the size L of the lattice, 
but still holding k = 0.1 fixed. We made L = 25, 50, 
100, 200, 300 and 400 for the application to the EW and 
KPZ equations. In the MH class, for which z — 4, we 
had to restrict our simulations to L — 20, 25, 30, 40, 50 
and 60, due to the large saturation times for larger sys- 
tems. The results are presented in figure where, in the 
left column, we show the roughness behavior for the sev- 
eral sizes L considered. In the right column we apply the 
Family- Vicsek scaling law @ , with the corresponding ex- 
pected value for a and z for each UC, in order to obtain 
the collapse of the various curves into a single one. The 
good collapses obtained, together with the results for (3, 
corroborate that our method indeed reproduces correctly 
each one of the three classes. 



t c ~ K K , 

t x - k- z * , (19) 

U sa t ~ K Q " . 

In the left of figureHHwe show the behavior of the rough- 
ness for different values of K, in the EW application of the 
method: for a system of size L = 250 we made k = 10 _1 , 
10~ 2 and 10 -3 , and averaged over 40 independent sam- 
ples. In the right of figure E3 we show the curve t c x n, 
where the power law fit provided z' K = 1.02(2). 

In figured we show the curves of saturation time and 
saturation roughness against k, in the EW application. 
Once again, the system size is L — 250 and the data is 
the average of 40 independent samples. As one can see, 
we found z K = 1.04(3) and a K = 0.509(2). 
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Figure 5: Log-log plots of the roughness for various system 
sizes L, in the three applications of the method (left column). 
As we apply the Family-Vicsek scaling law, using the expected 
exponents for each UC, good collapses are obtained (right col- 
umn). 
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Figure 7: Saturation time (left) and saturation roughness 
(right) as functions of the parameter k, for a system of size 
L = 250 and averaged over 40 samples. From the power law 
fits we obtained z K = 1.04(3) and a K = 0.509(2). 



k _1 plays the role of a characteristic time factor in the 
evolution of the system. Furthermore, z' K and z K must 
have the same value because otherwise, making k small 
enough, we would have either the uncorrelated regime 
taking place over the correlated one (for the case z' K > 
z K ), or the correlated behavior stretching over and over 
(for z' K < z K ), which cannot happen unless the system 
size is increased. In other words, the quantity t x — t c is 
supposed to be a function only of the system size L. 

The other result that we have obtained, a K = z' K /2, can 
be understood as follows. When t = t c , the roughness 
value is, say, u> — u> c and, considering that so far the 

1/2 

system is under the (3=1/2 regime, we have uj c = t c ■ 
It is also clear that uj c ~ uj sat , thus 

i *' 

UJ C ~ LO sat ti — LO sat => K~T — (t a « 




V. CONCLUSIONS AND PERSPECTIVES 



Figure 6: In the left panel, the temporal behavior of the rough- 
ness for L — 250 and k==10 _1 , 10~ 2 and 10~ 3 , averaged over 
40 samples, in the application to the EW class. The crossover 
between (5 = 1/2 (RD) and f) = 1/4 (EW) can be seen by 
comparing the curves with the dashed and dotted lines. In 
the right, the crossover time t c plotted against k, exhibiting a 
power law with exponent z' K = 1.02(2). 



For the other two classes, MH and KPZ, we have found 
very close values for the K-exponents as those obtained 
for the EW class. We conclude thus, that the crossover 
from random to correlated regime does not depend on 
the mechanism that generates correlations in the system. 
It is worthy to mention that for the KPZ class, when this 
crossover from random to correlated growth occurs, it is 
the Laplacian that dominates and the crossover should 
always be from random to linear correlated growth. 

The fact that we have found z K = z' K = 1 shows that 



We have introduced a new method, based on CA dy- 
namics, to study stochastic differential equations associ- 
ated to discrete deposition models. The method provides 
a new tool for obtaining the roughening exponents, which 
depends only on the discretization of the deterministic 
part of the equation, without having to actually solve it. 

We applied this method to study two linear equations 
(EW and MH equations) and a non-linear one (KPZ 
equation), in d = 1. The values obtained for the rough- 
ening exponents are in good agreement with prediction, 
showing that the method indeed reproduces each one of 
the three classes considered. In particular, for the non- 
linear case studied, a crossover from EW to KPZ class 
was obtained, for suitable values of parameter e, which 
controls the non-linearity strength. 

In addition, a crossover from RD to the considered cor- 
related class was obtained when we varied the parameter 
k. The crossover time, saturation time and saturation 
roughness were found to behave as power laws with k, 
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with numerical exponents z' K — 1.02(2), z K — 1.04(3) and 
a K = 0.509(2), respectively. These values have shown 
to be nearly the same, independently of the considered 
class. 

In further works, we proceed by applying the method 
to growth equations in which other terms appear, such 
as V 2 (Vft-) 2 and V • (V/i) 3 , which are the corrections up 
to fourth order to the V 2 h term in the EW equation |28l |. 



as well as verifying the validity of the method to growth 
processes in two-dimensional lattices, where discretiza- 
tion schemes are not as trivial as in one dimension. 
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